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Abstract 

The vast majority of sampling systems operate in a standard way: at each tick of a fixed- 
frequency master clock a digitizer reads out a voltage that corresponds to the value of some physical 
quantity and translates it into a bit pattern that is either transmitted, stored, or processed right 
away. Thus signal sampling at evenly spaced time intervals is the rule: however this is not always 
the case, and uneven sampling is sometimes unavoidable. 

While periodic or quasi-periodic uneven sampling of a deterministic signal can reasonably be 
expected to produce artifacts, it is much less obvious that the same happens with noise: here I 
show that this is indeed the case only for long-memory noise processes, i.e., power-law noises l/f a 
with a > 2. The resulting artifacts are usually a nuisance although they can be eliminated with a 
proper processing of the signal samples, but they could also be turned to advantage and used to 
encode information. 
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I. INTRODUCTION 



Nearly all digital signal-measuring equipment found in laboratories throughout the world, 
from the humble voltmeter up to powerful computerized data-logging systems and high- 
frequency digitizing scopes, operate with an internal master clock that sets the pace for an 
analog-to-digital converter that translates the electrical output of a transducer into a bit 
pattern. The master clock frequency is usually held as stable as possible and the sample 
intervals are fixed to a very high degree of precision. Correspondingly, most signal-analysis 
techniques are meant to be used on evenly spaced data: this is true for the Discrete Fourier 
Transform (DFT) and also for Autoregressive (AR) or Moving Average (MA) modeling of 
data [1]. However some data happen to be unevenly sampled: this is especially true for 
astronomers, who are seldom so lucky as to have an uninterrupted series of clear nights, 
and in general are bound to observe whatever comes from the sky, whenever it comes, and 
have to search for periodicities amid these scattered data. Indeed it was the astronomical 
community that developed the first effective spectral estimation techniques for unevenly 
sampled data M. 

n 

Uneven sampling has special properties: Beutler proved rigorously [3| that in general 
uneven sampling is not band-limited and later showed that a random, Poisson-distributed 

n n 

set of sampling times allows perfect signal reconstruction |^|. Earlier, Yen [5[ was able to 
derive modified forms of the Shannon reconstruction formula for different types of uneven 
sampling, which are however much more complex than the corresponding formula for even 
sampling. And indeed, randomly sampled signals are not easy to analyze and many standard 
methods must be abandoned, although in some cases one can restore regular sampling using 
reconstruction algorithms 

Since no sampling clock is quite perfect and is normally affected by noise 7] and by 
deterministic drifts (that may be periodic), all regular sampling should actually be regarded 
as quasi-regular sampling. Ignoring this may be dangerous, because it is clear that quasi- 
regular sampling of a deterministic signal (e.g. a sinusoidal signal) may introduce unwanted 
harmonics in the DFT analysis of the sampled signal, unless corrective measures are taken. 
But what happens if one samples pure noise? Can one still produce artifacts? If this were 
the case then an unrecognized quasi-regular sampling of a noisy background might become 
a problem since it would produce fake signals that could be mistaken for true. In this paper 
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I show that this is just what happens in some cases of colored noise, i.e., in the case of the 
long-memory noise processes l/f a with a > 2, while the correlation between samples for 
noises with a < 2 is insufficient to produce replicas of the low-frequency peak of the noise 
spectrum. A proof is given in section |H1 while section 11111 illustrates numerical results that 
confirm the theoretical analysis, and section HVI discusses some implications of these findings. 



II. PROOF BASED ON A DFT MODEL OF THE NOISE PROCESS 

It is well known that a signal sampled N times in the time interval (0, T) can be modeled 
by a sum of N exponentials and that this is equivalent to a DFT: in other words the DFT 



is a physical model of the signal [l] and we can write: 

1 ^— I „ / 2irinn, \ 



N ^ V N 

where f n denotes the n-th sample, and the fit coefficients correspond to the DFT and 
can be calculated from the formula 

^ 1 I „ ( 2irink\ ,„ s 

v 71=0 x ' 

With the usual regular sampling intervals At, the n-th sampling time is t n = nAt, the total 
sampling time is T = NAt, and equation (0) can be rewritten as follows: 

fn = -t= Fk exp ( 27[ik Y ) ^> 



k=0 

If sampling is not quite regular, the sampling times t n are replaced by t n + At n , and equation 
(JHJ) becomes 

, 1 c (2-Kink .k At n \ 

fn = — f= > F k exp — — h 2m— — — (4) 

VN _ \ n N At J w 

Now we assume a periodic sampling pattern with a period equal to M clock ticks, so that 

iV = mM and M < iV; we also assume at first that m is an integer (numerical simulations 

show that this requirement can be relaxed and m can be real) and we expand the relative 

timing shift At n / At clS db Fourier sum: 

At n 1 ^r-^ 1 / 2ninml \ 



1 r—s / 2ninml \ 

^Z^^exp — — (5) 



and in addition, we assume the relative timing shifts to be very small, i.e., At n / At <C 1. Next 
we notice that k/N < 1, and then, using equations (j3J) and (0) and after a few cumbersome 
but straightforward passages, we can approximate the observed DFT with the following 
formula: 

2 ' M ~ 1 

F ^ Fk + ~n~7W ^2 {k ~ ml ^ lFk ~ ml ( 6 ) 

1\ V M l=Q 

moreover if we make the rather weak assumption that the phase of the noise DFT in different 
frequency bins is uncorrelated so that (F k Fi) = if k ^ I where () is the usual ensemble 
average, then the DFT © gives the following spectrum: 
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l,l'=0 
o M-l 



1=0 

(the hypothesis of phase independence is quite common, because it is essential for noise 
generators like that of Timmer and Konig Q], and is supported by the numerical results 
reported in Q] ) - 

If the noise is white, i.e., the spectral density is flat, or if it is a 1/ f a noise with a spectral 
index a < 2, we see from eq. (J7J) that the periodic uneven sampling amounts to the addition 
of a (small) non-flat background. In fact a 1//* noise has a discrete spectrum S k ~ C/k a , 
therefore the observed spectrum (J7J) becomes 

. 2^ M-l 

s '^ Sk + w^H^ k - ml ) 2 ~ a \^ 2 ( 8 ) 

1=0 

However, if the spectral index a is greater than 2, then the l-th harmonic of the relative 
time shift At n /At produces a peak over the power-law background, which is just the low- 
frequency noise peak, shifted to the (m/)-th frequency bin. 
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III. NUMERICAL SIMULATION 



The analysis that leads to eq. J7J) assumes small relative timing shifts, but in this section I 
report numerical simulations carried out with the exact power-law noise generator described 
in P, U| that do support the analytical results also for large relative timing shifts 

The generator used in the simulation runs produces power-law noise from a superposition 
of random exponential pulses, and is exact in the sense that it produces a process that is 
theoretically guaranteed to yield a range-limited power-law spectrum between two extreme 
(angular) frequencies X min and X max . The generator takes correctly into account the cor- 
relation between samples in colored noises, and works also with unevenly spaced sampling 
times. 

In these simulations, time is in arbitrary units, and the average sampling interval is 
At = 1 arb. units; the choice of time units also sets the corresponding frequency units 
used for the relaxation rates X min and X max . Figure [T] shows a simulated signal obtained 
with the noise generator for a l// 3 noise; in this case the generator parameters are a = 3, 
^min = 0.0001, and X max = 1, i.e., the spectrum has a power-law region 1/ f 3 that spans the 
frequency interval 1.6- 10~ 5 < uj < 1.6- 10 _1 , and the pulse rate has been set at n = 10 pulses 
per unit time, so that the resulting noise signal is Gaussian to a very high degree |9j. In 
this case the sampling time has been sinusoidally modulated: At^/At = 1 + 0.2 sin(27rfc/4) 
(the period for uneven sampling is 4 samples), and there are in all 2 20 = 1048576 samples. 



Figure 121 shows the DFT spectrum of the signal o 
theoretical spectrum of the noise generator |ld 111 | 
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figure (0): a comparison with the exact 
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- which has a l// a power-law in the range A m j n < uj < X max - shows that on the whole 
the sampled noise process produced by the noise generator actually behaves as predicted by 
theory (9, 10, n|, except for a small peak at the frequency of the sampling time modulation. 



This small peak only shows up in this and in other runs (not shown here, but easily repro- 
ducible [12]) with a > 2, and this lends support to the proofs of the previous sections. The 
spectra shown in this and in the other figures have been partly detrended with a Hanning 
window (a general introduction to the need of the detrending step can be found, e.g., in jl^ ; 



see also the qualitative considerations in 




A closer look at the modulation peaks yields however a much more striking confirmation 
of the analytical results: in fact the theoretical spectral density (JOJ) of the noise generator 
has a 1/ f a power- law region for A m j n < uj < X max , while for uj < X m i n it has al// 2 behavior, 
and this means that from equation (J7|) we expect that the correction term due to uneven 
sampling is negligible just at the modulation frequency, while there should be two side- 
peaks whose exact shape depends on the low-frequency limit of the l/f a region, i.e., on 
\ m i n . And indeed this is just what happens in the simulations, as shown in figure El where 
part a. shows the region of the averaged spectrum in figure Eb close to the peak due to 
sampling time modulation. Figure EH> shows the expected behavior calculated from equation 
(J7J) and from the conditions used in the generation of the noise process and listed above (the 
correction for the incoherent gain of the Hanning window is also included). Finally, figure 
EH is the superposition of parts EH and EH ? where we see that the calculated shape closely 
matches the observed shape. 

Figure El shows the results of a similar calculation performed on the spectrum of a single 
realization of the noise process shown in figure EH: figure EH is the zoomed portion of the 
spectrum around the modulation peak. Figure Et> is a still closer zoom, and the arrow in 
the figure shows the position of the modulation frequency: notice that there is no peak 
just where one would naively expect to find one. The low frequency part of the spectrum 
in figure EH has been inserted in equation (|ZJ) to obtain the spectrum in figure EH, and we 
see that there is an almost perfect correspondence with the peak in figure Eh: this means 
that the correlation terms between different frequencies (the cross-terms in the derivation 
of equation (|7j)) are negligible even for a single realization of the noise process. 

I have also noted that there must be a dependence of the split-peak shape on the exact 
shape of the low- frequency part of the spectrum: figure EH shows the shape of the peak for a 
larger value of \ m in (A m m = 1CT 3 ). In this case the low-frequency part of the noise spectrum 
(EJ) has a wider 1/f 2 region, and the side-peaks must be correspondingly lower and further 
apart: indeed this is just what happens in figure EH- A further confirmation is provided by 
figure Eb, which shows the peak for a smaller value of X min (A min = 5 ■ 1CT 5 ): the side-peaks 
are much higher and also closer. 
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IV. DISCUSSION 



While most of the observed power-law (1/ f a ) noises have spectral indexes < a < 2, with 
an apparent clustering around a — 1, red noises, i.e., noises with spectral indexes a > 2, also 
show up in several unrelated systems 16, 17] like the water level of the Nile river, economics, 

n 

orchid population size |18|| and local temperature fluctuations and affect precise timekeeping 
jl^ and our ability to predict environmental and animal population variables!^. Noises 
with a > 2 also appear in the energy level fluctuations of quantum systems |2l|, |22j and in 
timing noise in pulsars j^J. Because of their extreme peaking behavior at low frequencies 
these noises are also called "black" jl7j, and they display marked persistence properties 1161 
that may lead to the mistaken identification of underlying trends in experimental data [24]. 
From the results reported in this paper it follows that these noises pose yet another potential 
danger to experiments that use uneven sampling, because their long-memory properties give 
rise to artifacts in the DFT spectra. 

However equation (|ZJ) also shows that the single peak - in the case of true l/f a noise 
- or the side-peaks - in the case of range-limited l/f a noise - can be modulated both in 
amplitude and in frequency by modulating either the noise spectrum or the relative timing 
shift amplitude, or the repetition index m: this means that uneven sampling of colored noise 
could be utilized to encode information, and since an encoding noise appears at first sight 
just ordinary noise, this could be used to implement a secure communication channel (there 
is a very rich literature on this topic, but here I give only a reference to a classic book 
25I ] and to a recent paper (2^). A simple example of the kind of modulation that can be 
achieved can be gleaned from figures E3 ID and EJ if one uses the noise generator J9, 10, 11 1, it 
is possible to modulate the shape of the low-frequency part of the spectrum with a proper 
change of X m in and in this way one modulates in turn both the amplitude and the position, 
i.e., the frequency, of the side peaks. 
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FIG. 1: This figure show a noise signal produced with the generator described in 
The parameters in this specific run are a = 3, \ m in = 0.0001, and A max = 1, i.e., the spectrum 
has a power-law region 1/f 3 that spans the angular frequency interval \ m in < w < A max . The 
generator produces power-law noise from a superposition of random exponential pulses, and in 
this run the pulse rate has been set at n = 10 pulses per unit time, so that the resulting noise 
signal is Gaussian to a very high degree [{J. Time is in arbitrary units, and the average sampling 
interval is At = 1 (arb. units); the choice of time units also sets the corresponding frequency units 
used for the relaxation rates A m , n and A max . The sampling time has been sinusoidally modulated: 
Atfc/Ai = l + 0.2sin(2-7r£;/4), and the signal has been sampled 2 20 = 1048576 times. Part a. shows 
the whole signal generated in this run (time does not start from zero, because at the beginning 
some samples are used for the generator initialization and are discarded); part b. shows the initial 
32768 valid samples. Notice that even though the sampling time modulation is rather large (20%), 
it is quite invisible in the zoomed figure. 
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FIG. 2: a. Spectrum of the signal shown in figure Q vs. the scaled angular frequency u/vg- The 
arrow marks the peak from sampling time modulation at the expected position uj = wtv/2, where 
ujn = is the (angular) Nyquist frequency and v$ is the sampling frequency, b. Spectrum aver- 
aged over 16 realizations of the same noise process: the solid line shows the expected (theoretical) 



behavior of the noise spectrum 
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111 ] , corrected for the incoherent gain of the Hanning window 



that has been used for trend removal. The upward bend at high frequency in the spectrum b is 
due to aliasing which becomes apparent after averaging, while the upward bend at low frequency 
is due to the uncorrected DC component which cannot be removed by windowing. 
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FIG. 3: a.) Averaged spectrum of figure zoomed and centered on the position of the small 
peak (i.e. u)/u$ = 7r /2), with linear scales on both axes, b.) Expected behavior calculated from 
equation (J7J and from the conditions used in the generation of the signal in figure ^ (the correction 
for the incoherent gain of the Hanning window is also included), c.) Superposition of parts a and 
b: we see that the calculated shape closely matches the observed shape b. 
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FIG. 4: a.) Spectrum of figure [2K zoomed and centered on the position of the small peak (i.e. 
lu/vs = tt/2), with linear scales on both axes, b.) An even closer view of the peak shows that it 
is actually a split peak: the arrow marks the position of the modulation frequency, c.) This part 
shows what one obtains if one uses equation (JJJ) and the low-frequency part of the spectrum of 
figure^ to calculate the expected shape of the peak. The calculation includes the correction for 
the incoherent gain of the Hanning window that has been used for trend removal. We see that the 
calculated shape c. matches very well the observed shape b. 
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FIG. 5: This figure illustrates how the peak splitting changes for different noise shapes: a.) In 
this case the noise generation parameters are the same as for the signal in figure ^ except for 
the lowest relaxation rate \ m i n = 0.001, and the figure shows the averaged spectrum zoomed and 
centered on the position of the small peak. The softer low-frequency behavior in noise spectrum 
produces side peaks that are smaller and further apart than those shown in figure The solid line 
shows the expected behavior, calculated as in figureOJa b.) A smaller value of the lowest relaxation 
rate, \ m i n = 0.00005 yields instead much closer and higher side peaks. Once again we see that the 
calculated shape (solid line) matches very well the observed shape. Because of the much narrower 
splitting, the scale in part b is expanded with respect to a to improve visibility. 
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